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The Gaussian state description of continuous variables is adapted to describe the quantum in- 
teraction between macroscopic atomic samples and continuous-wave light beams. The formalism 
is very efficient: a non- linear differential equation for the covariance matrix of the atomic system 
explicitly accounts for both the unitary evolution, the dissipation and noise due to the atom-light 
. . . interaction, and the back-action due to homodyne optical detection on the beam after its interac- 

' tion with the atoms. Applications to atomic spin squeezing and estimation of unknown classical 

I parameters are presented, and extensions beyond the Gaussian states are discussed. 
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I. INTRODUCTION 



Pulses of light, large atomic ensembles, and collections of more than, say, hundred trapped ions, are quantum systems 
where the behavior of various collective degrees of freedom is well described by quantities which have continuous 
I ' spectra, i.e., the systems may be described by collective effective position and momentum variables. The demonstration 
^ , of quantum control of these systems varies from studies of squeezing and entanglement, over storage and retrieval of 
\l optical information in gases to high precision probing of classical properties in atomic magnetometry, atomic clocks 
] and inertial sensors. Control is exercised via tunable interactions, by state reduction due to measurements on the 
^ ■ systems, and by feed-back schemes in connection with measurements. 

I [ In quantum optics, the quantum properties of a continuous beam of light arc normally described in the Heisenberg 
iy~j i; picture, where field operators are expressed (often in the frequency domain) in terms of incoming vacuum fields 
■ with standard correlation functions. This input-output formalism leads, e.g., to the noise spectrum of a squeezed 
, light beam[p. This approach accounts for the results one obtains if measurements are carried out directly on the 
^ ^' beam, but it has been technically very difficult to describe the situation where the light beam is made subject to 
^ I ^ ] interaction with another quantum system and is subsequently measured. The measurement record is stochastic, and a 
^ ■ real-time description of the measurement back-action on the probed quantum system is normally referred to quantum 
' trajectory or Monte Carlo wave function treatments in the Schrodinger picture, which are incompatible with the 
^ , frequency domain Heisenberg representation of the optical beam. 

^ '' In Fig.^ we display the interaction between light and atoms. A Gaussian state analysis was introduced recently to 
deal, in general terms, with the quantum properties of these systemSjand tools were developed to handle interactions 
and measurements which preserve the Gaussian state character 0, 13 . As we shall illustrate below, the Gaussian 
^ . ^ description is useful because (i) it handles the interaction between atoms and a quantized continuous-wave (cw) beam 
' of light and (ii) it allows a description of measurement induced back-action in real time. This description thus provides 
a useful approach to a long standing problem in quantum optics, and it presents a theoretical treatment of physical 
systems and interactions of high current interest. The description is restricted to Gaussian states. A cw laser beam 
described by a coherent state and squeezed and quantum correlated optical beams created by down conversion are 
Gaussian in the field canonical variables and hence readily incorporated in our treatment. Turning now to the atoms, 
our approach does not describe the interaction with a single ion or atom, but a very accurate mapping exists between 
macroscopically spin-polarized atomic samples and a single harmonic oscillator. 

In this work, we describe the practical application of the Gaussian state formalism to continuous variable systems, 
allowing full account of back-action due to measurement, noise, losses and inhomogeneities of the systems. The for- 
malism is illustrated by a discussion of explicit examples concerning spin squeezing, magnetometry and entanglement. 
A whole tool-box can be created, describing the effect of frequency filters, finite band-width sources and detectors, 
finite efficiency detection, and dark counts, simply by adding extra reservoir modes. In practice, the Gaussian state 
for a system of n quantum harmonic oscillators, representing a number of optical beams and atomic components, is 
described by 2n mean values for the quadrature components and by a 2n x 2n covariance matrix. While the evolution 
during measurements of mean values is stochastic, the covariance matrix is propagated in time in a deterministic way 
(see Sec. 2). It is a remarkable advantage of the Gaussian state description that extra physical systems and reservoir 
modes can be included at only little expense (two extra rows and columns in the covariance matrix per mode). In the 
last paragraphs of this work, it will be discussed how to develop a theory for continuous variable systems where the 
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FIG. 1: Atom-light interaction. In the figure, we display a cloud of atoms described by collective continuous variables iatjPat, 
interacting with (a) a pulse of light and (b)-(c) a continuous wave of light. In (a) the assumption of an accurate description of 
the light pulse in terms of just a single mode (aiphiPph) leads to a simple and natural description in the time domain of both the 
interaction of the pulse with the atoms as well as the detection process. In quantum optics, the continuous beam of light in (b) 
is normally described in the frequency domain, say, by canonical operators Xph,^^, j5ph,t^. The interaction with the atoms and the 
measurement process, on the other hand, is more readily described in the time domain, and as discussed in detail in the text, 
it is technically difficult to pass from the frequency domain for the light operators to the time domain for the description of 
interaction and measurement process. To circumvent this problem, we introduce an effective description of the integral system 
in the time domain as indicated in (c). Here the beam is divided into segments of duration t and length L = ct each of which 
is assumed to be short enough to be accurately described by a single mode Xph,i,Pph,i, and the interaction with the atoms and 
the measurement is described by a succession of interactions with the individual beam segments. 



Gaussian description breaks down, either because of the interactions involved, because of the measurement schemes, 
or because of couphng of a small discrete system to collective continuous degrees of freedom. 



II. TIME EVOLUTION OF GAUSSIAN STATES, GENERAL THEORY 



In this section, we introduce the Gaussian description in a general setting using existing results [2, lilj]. Sections 2.1 
and 2.2 deal with the evolution of continuous variable systems due to a bilinear Hamiltonian and linear losses. This 
evolution can be solved by an affine transformation in time of the canonical operators, and all system properties are 
given by their mean values and their covariance matrix for which an exact treatment is provided. Section 2.3 deals 
with the effect of measurements on the system. The update of the system state vector or density operator conditioned 
on a measurement outcome is non-trivial in the most general case, but as we shall see, Gaussian states transform into 
other Gaussian states in a well described manner under homodyne detection on part of the system, and in this case 
the mean values and the covariance matrix still provide all properties of the system. 



A. Time evolution due to a bilinear Hamiltonian 



Let y = {xi,pi,X2,P2, ■ ■ ■ ,Xn,'Pn)'^ denote the column vector of 2n variables with canonical commutators = 
iSij, and let H = H{xi,pi,X2,P2, ■ • ■ , Xn,Pn) denote the Hamiltonian of the system. We shall assume that H is bilinear 
in the canonical variables. Heisenberg's equations of motion during time r are then solved by a linear transformation 
of the operators by the matrix 



y{t + T)=^ry{t)- 



(1) 
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The same transformation applies to the vector of mean values m = (y), m(t + r) = S^m(t). From Eq. |^ and the 
definition of the covariance matrix 7^ = 2Re ((y^ — {yi)){yj — (j/j))), we directly verify that 7 transforms as 

7(t + r) = S,7(0S^ (2) 

under the interaction. 



B. Time evolution due to dissipation and noise 

In the absence of dissipation Eq. (j^J determines the evolution of the covariance matrix. In realistic situations, 
however, there will be sources of dissipation and noise. Dissipation leads to a reduction in the mean values of the 
canonical variables, and as is known from the quantum theory of damping and the fluctuation-dissipation theorem of 
statistical mechanics, such a reduction must be accompanied by fluctuations. In the quantum domain we must, e.g., 
fulfill the Heisenberg uncertainty relations, also when the mean values are reduced. The generalization of Eq. ||2J) to 
the noisy case reads for small r 

7(t + r) =L,S,7(t)S^L,+N,, (3) 

where L,- describes the reduction of the mean values, m(t + t) = Lr^r^{t), and where Nt- is the associated noise. 
In examples below, we shall give explicit forms of these matrices. 

If the state of the system is initially a Gaussian state, i.e., its Wigner function for the canonical variables is a 
Gaussian function, the evolution due to a bilinear Hamiltonian preserves the Gaussian character. The same is true 
for linear damping of an optical field mode, and as validated by a calculation and more detailed discussion^ it also 
holds to an excellent approximation for atomic decay models. 



C. Time evolution due to a homodyne measurement event 



The above arguments were based on the Heisenberg picture evolution of the canonical operators, but the evolution 
due to measurements is more conveniently described as state reduction in a Schrodinger picture representation of 
the system state vector or density operator. A general representation of the state, pure or mixed, of a collection 
of harmonic oscillators is provided by the Wigner function W{£,) with ^ = {£,1, ■■■£,2n) S R^". This function is 
connected with the density matrix in position or momentum representations by a Fourier-transformation, and it 
provides a good intuitive picture of the phase space distribution of the system. In fact, the expectation value of any 
symmetrically ordered function Fsym{S:i,Pi, ■■■Xn,Pn) (-Fsym is the average of all the ways of ordering the operators 
defining F{xi,pi, ...a;„,p„)), is given by the pseudo-classical expression: 

(i^sym(xi,Pl,...£„,p„)) = / d^"^ WiOFiO- (4) 



We recall that we aim at a description of the state of an atomic sample subject to interaction with an optical beam 
which is being probed after the interaction. We hence address what happens to the quantum state of the remaining 
system when one of the sub-systems (with a conjugate pair of observables Xn,Pn) is subject to a measurement. 
Examples of measurements are positive operator valued measures with coherent state outcomes, homodyne detection 
which projects the measured sub-system onto a position or momentum eigenstate (equivalent to the limit of a strongly 
quadrature squeezed state), and number state detection. Such measurements project the (xmpn) sub-system onto 
a particular state which we can also describe by a Wigner function Wmceis{^2n-i, ^2n)- The state of the remaining 
system conditioned on the outcome leading to this particular state is 

^VcondUl, • • • ,42n-2j = ^ , (o) 

with Pmoas = J d^l... d^2nyV{^l, . . . , 6" ) Wmoas (6n-l , C2«) • 

Now turning to the Gaussian states, a series of simplifications occur. For example, the Wigner functions 
W'mcas(C2n-ii C2n) for cohcrcut and squeezed states are Gaussian functions of the variables. This implies, that if 
the initial Wigner function is a Gaussian function of the variables, this property is maintained by the homodyne 
detection process. Generally, the Wigner function for a Gaussian state is fully parameterized by the mean values m 
and the covariance matrix 7: 

WGauss(e) = exp (-(e - m)^7-i(C - m)) . (6) 
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As the Gaussian character is also maintained by the bilinear Hamiltonian and the linear decay processes, we conclude 
that to describe the time evolution of a system which starts in a Gaussian state, it suffices to provide the time 
dependent m and 7. 

Since part of the system is being measured upon, and hence disappears from our quantum state, cf. Eq.jSJ, it 
makes sense to write the covariance matrix in the form 

where the (2?! — 2) x (2n — 2) sub-matrix A is the covariance matrix for the variables yi — {xi,pi, . . . , Xn-i,Pn-i)'^ 
which are not subject to measurement, B is the 2x2 covariance matrix for the sub-system subject to measurement 
y2 = (irijPn)"^! and C is the 2 x (2n — 2) correlation matrix between the elements of yi and 2/2- According to the 
above expressions Q-©, a measurement of x„ transforms A as 2, ^3, 3 

Ak^ A' = A-C(7rB7r)^C^, (8) 

where tt = diag(l,0), and where ( )~ denotes the Moore-Penrose pseudoinverse: (ttBtt)" = diag(i3(l, 0). If we 
associate with the precise measurement of Xn an infinite variance of pn and hence a total loss of correlations between 
Pn and the other observables, this result is equivalent with the Bayesian update of a classical Gaussian probability 
distribution We recognize the Moore-Penrose pseudoinverse as the normal inverse of the corresponding covariance 
matrix, (ttBtt)^ =diag(i?(l, 1), oo)""'^. 

Unlike the covariance matrix update, which is independent of the value measured, the vector m = (y) of expectation 
values will change in a stochastic manner depending on the actual outcome of the measurement. The outcome of the 
measurement of Xn is random, and the measurement changes the expectation value of all other observables due to 
the correlations represented by the covariance matrix. Let x denote the difference between the measurement outcome 
and the expectation value of in, i.e., a Gaussian random variable with mean value zero and variance given by half of 
the appropriate covariance matrix element B{1, 1). It follows again from Eqs. (O-®, (and from the corresponding 
classical theory of multi-variate Gaussian distributions,) that the change of mi = (yi) due to the measurement is 
given by: 

mi m; = mi + C^(7rB7r)-(x,-)^, (9) 

where we use that {■kBtt)~ = diag(_B(l, 1)"^, 0), and hence the second entrance in the vector (x, ■) need not be 
specified. 



D. Time evolution due to continuous homodyne measurements 



In the continuous interaction between a cw light beam and a cloud of atoms one faces a situation where a single 
system (the atoms) is continuously indirectly monitored, e.g., by a homodyne detection of the light field. This raises 
the problems, mentioned in the introduction, of dealing simultaneously with a continuous beam and measurement 
induced back-action (see also Fig.^). We have recently solved this prolslem for Gaussian states0 by quantizing the 
light beam in short segments of duration t and corresponding length L — ct. These beam segments are chosen so 
short that the field in a single segment can be treated as a single mode and such that the state of the atoms interacting 
with the field does not change appreciably during time r. The evolution of the atomic system with the entire beam 
of light is obtained by sequential interaction with subsequent light segments. The generic multi-mode character of 
the cw beam of light is treated in the Schrodinger picture in time domain rather than in the Heisenberg picture in 
frequency domain (cf. Fig. 1(c)). 

The simplest example of continuous light-atom interaction is the one of a coherent monochromatic beam of light, 
corresponding to a product state of coherent states in each segment along the beam axis. In this case, the problem 
simplifies significantly because all segments are in the same trivial state prior to the interaction with the atoms. The 
segments need not be included formally in the update of the covariance matrix until it is their turn to interact with 
the atoms. Segments which have already interacted with the atoms may be detected instantly after the interaction, 
and in practice they are if the detector is placed within meters from the interaction volume. The detected segments 
then disappear from the formal description of the system. Prior to the interaction with the beam, we thus consider 
only the atomic covariance matrix A, and in the absence of any correlation with the incident beam segment, the 
block-off-diagonal matrices in Eq.(|7I) vanish 



C = 02x(2n-2) 



(10) 
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while the field state of the incident segment is characterized by the normal noise properties of the coherent state 

B=l2x2. (11) 

The full covariance matrix is now propagated according to Eq. iPJ , and the matrix changes to describe the state of 
the atoms and the optical segment after interaction. To describe the effect on the atoms of the measurement on the 
field segment, we apply the measurement update formula (jSl for the atomic part, and since the field segment has been 
observed and reduced to classical information, we are ready to turn to the interaction with the next light segment, 
which conveniently fits into the covariance matrix in the same locations as the previous segment according to 
Eqs. (|10l) - (|ll|l . This evolution is repeated to describe in real time the interaction with a beam for any extended 
period of time, and the expectation value and our uncertainty about any variable of the system at the end of the 
interaction is readily found from the appropriate entrances in the vector m and the matrix 7. 

In the limit of small t the changes in 7 and m expressed by the update formulae ©, and are infinitesimally 

small. In this, suitably defined, continuous limit, the update formulae translate into differential equations. After 
application of Eq. Q, the sub-matrix C depends linearly on the elements of A and as shown in Eq. H21() and the 
ensuing discussion below, its elements are proportional to ^/t. B is essentially unchanged for short t, and A changes 
linearly with r. In the limit of infinitesimally small time increments, the update formula may therefore be written as 
a closed non-linear equation of motion for A: 

A = lim ^ ^ = G - DA - AE - AFA, (12) 

with suitably defined matrices G, D, E, F. This equation is an example of a so-called matrix Ricatti equation(8i|, and 
by the decomposition A = WU^^, it can be rewritten in terms of two coupled linear equations W = — DW -I- GU, 
and U = FW + EU. Below, we shall see examples of analytical solutions to the problem based on these equations. 

III. APPLICATION OF THE GAUSSIAN FORMALISM TO ATOM-LIGHT INTERACTION 

The Gaussian formalism can be applied to describe the interaction between atomic samples and optical beams. In 
our examples, we consider optical Faraday rotation, which probes the collective spin ground state of a gas of atoms. 
To introduce the transition to an effective Hamiltonian expressed in terms of canonical variables, we discuss in some 
detail the interaction of an atomic ensemble with a pulse, or segment, of light. 

A. Stokes vector and canonical conjugate variables for light 

To make the discussion simple, at first only a single atomic sample and a single pulse or segment of a light beam will 
be considered. In Faraday rotation experiments, one uses light, which is linearly polarized along the, say, x-axis. The 
interesting quantum degree of freedom of the light pulse is not the field amplitude itself, but the intensity difference 
between the linearly polarized components along 45 and 135 degree directions in the xy plane, and between the two 
circularly polarized components with respect to the z-axis. These components are equally populated on average, 
but as every single cc-polarized photon can be expanded as a superposition of single photon states of either pair of 
polarizations, their populations will fiuctuate according to a binomial distribution. For a pulse with a definite number 
-^ph of photons, one may represent these populations conveniently by the components of the Stokes vector, where the 
cc, y and z-components represent the populations difference of x and y polarizations, 45 and 135 degree polarizations 
and tT+ and ct^ -polarizations, respectively, i.e., 

= ^ {ala,^ - at ay) = - ^ {h\a^ + a^a+) , (13) 
^^^ \ ("J^fly + 4^^) = (a^a_ - a}_a^ , (14) 
'S'^ = ^ {aiay - ala^) = | (^a^a+ - ala-) . (15) 
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Since the light is assumed to be Hnearly polarized along the x axis, Sx niay be treated classically and from Eq. H13|l . 
Sx/f^ = Sx/fi = A^ph/2. The Stokes vector components obey the commutator relations of a fictitious spin, and the 
variance of the binomial distributions are in precise correspondence with the quantum mechanical uncertainty on Sy 
and Sz^ achieving the Heisenberg hmit Ya,Y{Sy)Ya.Y{Sz) — \{hSx)\'^ 1^- 

We assume that Sx remains large and essentially unchanged during the interaction with the atomic gas, and we 
can then introduce the effective position and momentum operators 

(iph,Pph) = (Sy/^/W^\,S,l./W^\) , (16) 

which fulfill the standard commutator relation [iphiPph] = i and resulting uncertainty relation. These are the canonical 
conjugate variables that we wish to describe by the formalism outlined in the previous section. The initial binomial 
distributions of 5*^, Sz approach Gaussian distributions in the limit of large photon numbers. Moreover, the fact that 
the uncertainty relation is minimized in the initial state implies that this state is a Gaussian state, i.e., the Wigner 
function for the field is a Gaussian function l9||. 



B. Atom-light interaction 

The physical system of interest consists of one or more macroscopic ensembles of trapped atoms interacting off- 
resonantly with one or more laser beams. We consider the usual electric dipole interaction between the atoms and the 
quantized field. First, the off-resonant coupling of the atoms with the light field is expanded in transition operators 
between the ground {\FM)) and excited {\F' M')) hyperfine states (several excited states with different F' may be 
coupled to the ground state) . Then, the atomic coherences pertaining to the excited states are expressed by the light 
fields and ground state coherences by adiabatic elimination using Heisenberg's equations of motion for the slowly 
varying operators. This procedure generally allows us to derive a dispersive effective Hamiltonian|lfl lllj|. which for 
the Na^t atoms reads [T^ 

Wat F 

-ffint.r [(c+,Af(A)a^a+ +c_,M(A)aU_) \FM)j{FM\ 

j=l M=-F 

+bM{^) (^ala+\FM + l)j{FM - 1\ + a^a_|FAf - l)j{FM + 1|)], (17) 

where field creation and annihilation operators for (t+ and a^-polarized photons have been introduced. The first two 
terms describe the ac Stark shift of the ground state \FM) caused by the coupling to the excited \F'M ± 1) states by 
the two field components. The coupling coefficients are given by c±^m(A) = -''2h'^p,{gpJ^.l.p,J^,)'^/AFr where Ap> is 
the detuning of the laser frequency from the upper level, and where the coupling constants gpM-F'M' electric 

dipole coupling matrix elements, dpM-F'M' ~ ^/^^^^oT^^^^^o^^d^ m-f' m' ■ These matrix elements contain the 'electric 
field per photon' for a plane wave field with transverse area A and length cr, and they involve the spherical tensor 
components of the dipole operator d — — er of the electron, dpj^.p,j^,j, = {FM\d±\F' M') and d+ = —{dx + idy)/\^, 

(i_ = {dx — idy)/^/2. The terms in Eq. H17() proportional to 6m(A) ~ '~'^f^J2F' 9fm-i-f'm9fm+i-f'm/^P' describe 
AM — ±2 Raman transitions involving absorption and stimulated emission of a pair of photons with different 
polarization. 

1. Spin 1/2-case 

For much of the discussion in the rest of this work, we shall restrict ourselves to the case of atoms with only one 
ground and one excited level which both have total angular momenta F — F' — 1/2. The above Hamiltonian simplifies 
in this case, and noting further that the dipole matrix elements are related to the total spontaneous decay rate F of 
the upper state, c±.zpi/2 = — 3?iF(t/(2tAA), with a = A^/(27r) the resonant photon absorption cross section, Ea. H17|) 
then reduces to 




where the F — F' = 1/2 index has been suppressed. 



7 



The atomic ensemble is initially prepared with all iVat atoms in a superposition (| — 1/2) + |l/2))/-\/2 of the 
two ground states with respect to the quantization axis z, i.e., the total state of the atoms is initially given by 
((I — 1/2) + |l/2))/V2) In this state, the system of two-level atoms is described by a collective spin, «/ = f J2j 
with a component along the x-direction which attains the macroscopic value (J^;) — ^-^at/2, and with a component 
along the z-axis, J^, which represents the population difference between the | ± 1/2) states. Similarly, we may use 
Eq. lfT5|l and represent the operators of the photon field in terms of the collective Stokes vector operator, Sz- The 
Hamiltonian can, hence, be rewritten in terms of the collective spin variables for photons and atoms 

Hint,T = T—rSzJz, (19) 

tAA 

where an overall energy-shift proportional to the number of photons in the pulse segment has been neglected. 
As for the photons it is convenient to introduce effective atomic position and momentum coordinates 

(iat,Pat) - (^Jy/\/\{flX)\, Jz/^\{hX)\j , (20) 

for which the initial state is a minimum uncertainty Gaussian state. The last step of this analysis is then to rewrite 
Eq. H19(l in terms of canonical conjugate variables, 

-ffint.r = ^KrPatPph, (21) 

where 

= -:^V&V^^^- (22) 

The Hamiltonian correlates the atoms and the light fields and is bilinear in the canonical variables. Hence the 
theoretical formalism of Sec. 2 applies. The coupling constant is small for realistic parameters, and a coarse 
grained description, where the atoms interact with one segment of light after the other, will be perfectly valid even 
for the macroscopic number of photons iVph in each segment required by our Gaussian treatment. Note that {S^) is 
proportional to the number of photons in the beam segment, i.e., to r, and it follows that Hiat,TT is proportional to 
y/r yielding a well-defined differential limit in Eq. p2(l . 

We have emphasized the convenience of using Gaussian states, because their Schrodinger picture representation is 
very efficient and compact. Now, given that every segment of the optical beam becomes correlated with the atomic 
sample, as a function of time, the joint state of the atom and field has to be specified by a larger and larger number 
of mean values and second order moments. If no further interactions take place between the atoms and the light after 
the interaction, there is no need to keep track of the state of the total system. In practice, either the transmitted 
light may simply disappear or it may be registered in a detection process. In the former case, the relevant description 
of the remaining system is obtained by a partial trace over the field state, which produces a new Gaussian state of 
the atoms, which is simply given by removing the photonic lines and columns of the covariance matrix immediately 
after the interaction update The measurement of the small Faraday rotation of the linearly polarized probe is 
done by a measurement of the intensity difference between the 45 and 135 degree polarization components, i.e., by a 
measurement of the Sy oc Xph observable, which is precisely the "homodyne" measurement described in section FU PI 
The atomic state is thus described by the corresponding update formula of Eq. H12|l . 



IV. SPIN SQUEEZING IN THE GAUSSIAN DESCRIPTION 

With spin squeezed atomic ensembles, i.e., samples where the variance of one of the angular momentum (spin) 
components is reduced compared with the coherent state value, one has the possibility to measure certain atomic 
and/or classical parameters beyond the precision set by the standard quantum noise. 

The theory of squeezing of the collective atomic spin variable was dealt with in a series of papers [13, [Hill, and 
extended to include investigations of quantum non-demolition feedback schemes|l4l Il5j . and inhomogeneous light- 
atom coupling[l2, [I3|- In a series of related works[3 [l8|, (l^ |20|, spin-squeezing of continuous variable quantum 
systems has been investigated in the approximation where the atomic and photonic degrees of freedom are described 
by a Gaussian state. 

We are interested in the case, where the polarization rotation of the light field is registered, i.e., the observable 
iph is measured. The effect of measuring one of the components in a multi-variable Gaussian state is effectively to 
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produce a new Gaussian state of the remaining variables as discussed in detail in Sec. The column vector of the 
variables for the gas and the photon field reads y — (iatjPat, XphiPph)-^ and the S'-matrix in Eq. is 

/I Q Kr 
„ _ I 1 

^ ^ 1 

\ 1 




A. Dissipation and noise 

In the probing process there is a small probability that the excited atomic levels which were adiabatically eliminated 
from the interaction Hamiltonian of Eq. 1)21(1 will be populated. If this happens, the subsequent decay to one of the 

two Mz — ±1/2 ground states occurs with the rate 77 = ^ ) ' '^ti^^^ ^ is the photon flux and where the 
remaining parameters were defined in Sec. IIII B II The consequence of the decay is a loss of spin polarization since a 
detection of the fluorescence photons in principle can tell to which ground state the atom decayed. If every atom has 
a probability rjr — rjr to decay in time t with equal probability into the two ground states, the collective mean spin 
vector is reduced by the corresponding factor (J) (•/)(! — ?7t)- 




t [ms] 

FIG. 2: Uncertainty of pat as function of time during optical Faraday rotation experiment. The effective coupling is = 
1.83 X 10S-\ The lower curve is without inclusion of atomic decay, and the upper curve includes atomic decay with a rate 
77 = 1.7577 s~^ and photon absorption with a probability e = 0.028. These values correspond, for example, to a 2 mm'^ 
interaction area, 2 x 10^^ atoms, 5 x 10^'* photons s~^, 10 GHz detuning, and 852nm light, appropriate for the ^^^Cs(6S'i/2(-F — 
4) — 6Pi/2(-F = 5)) transition. Factors of order unity related to the coupling matrix elements among different states of the 
actual Zeeman substructure of Cs are omitted. 

Simultaneously, every photon on its way through the atomic gas has a probability for being absorbed|l9| e = 
TyJ+/\-i ) (see Sec. IIII B"T1 for definition of parameters). The effect of these noise contributions were discussed 
in detail elsewhere [sl l20j . and the result for the reduction and noise matrices of the update formula of Eq. reads 
= diag(Vl - Vr, VI - Vr, VI e), and = dia.g{-^^^ r]r, -^^Vr, 2(s!'(t)) '^' 2(5!^) ^^ ^^"^ ''^''^ ^ ^■ 

The factor hNut/ {Jx{t)) initially attains the value 2, and increases by the factor (1 — r/r)^^ in each time step r. The 
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factor hNpii/ (2{Sx{t))) is initially unity, and is approximately constant in time since the light field is continuously 
renewed by new segments of the light beam interacting with the atoms. 

We note that when the classical x-component of the atomic spin is reduced this leads to a reduction with time of 
the coupling strength — rjr (see Eq. (|22f) '). 



B. Solution of Ricatti equation 



We now have explicit forms for the matrices needed for our update of the Gaussian states. In the Gaussian 
description, the problem of spin squeezing may be solved either by the discrete update formulae or analytically from 
the matrix Ricatti equation. In the latter case, we note that the covariance matrix after n iterations in the noise-less 
case is 



7n 



2Var(xat) 

2Var(pat) 

10 

1 , 



We then apply the S'-matrix from Eq. (|23|l and find 



(24) 



2Var(xat)+K^ K-r 

2Var(p,t) 2K^Var(33^t) 

2K^Var(33at) 1+2k^ VarCjS^t ) 











(25) 




(2Var(35at))" 



2 ] , insert into Eq. ||SJ|, take the 



From this matrix, we determine, to lowest order in t, C(7rB7r) C = n 
continuous limit and use — k^tIt- This procedure leads to the following differential equation for the variance of 
Pat(oc Jz)'- ^Var(pat) — — 2^2 (Var(pat))^, which is readily solved by separating the variables 



Var(pat) 



2kH + l/Var(pat,o) ' 



(26) 



where Var(pat,o) = 1/2 is the variance of the initial minimum uncertainty state. Note that the solution to the variance 
of the conjugate atomic variable is Var(a;at) = K^tj^ + Var(a;at,o) with Var(i;at,o) = 1/2. Hence, while Pat is squeezed, 
Xat is antisqueezed to maintain the equal sign in Heisenberg's uncertainty relation. 

When dissipation and noise is included the problem may still be solved analvticallv j2n|| . The expressions for the 
variances are quite complicated and will not be given here. Figure |21 shows the spin squeezing as a function of probing 
time. When atomic decay is not included, the uncertainty in Pat is a monotonically decreasing function with time. 
When decay and noise is included, a minimum at t„iin is reached whereafter the degree of squeezing starts to decrease. 
On the time scale of the figure, which is chosen to reflect realistic experimental time scales, the increase in Var(pat) 
is hardly visible. 



C. Inhomogeneous coupling 

One of the virtues of the Gaussian description of spin squeezing is that it is straightforwardly generalized to handle 
situations which are hard to approach by standard means. For example, a variation in the intensity of the light 
beam across the atomic sample and a large photon absorption probability both lead to an inhomogeneous atom-light 
coupling^, |23. To treat such a case, the atomic gas is divided into n slices each with local light-atom coupling 
strength k.^. The 2n -f 2 dimensional vector of gaussian variables describing the 2n collective canonical position and 
momentum variables for the atoms, and the two collective position and momentum variables for the photon field then 
reads y = (iat.iiPat.i, • ■ • , £at,n,Pat,n, iph,Pph)^, and the generalization of Eq. to this case is 

ffint.r = fi- K^.iPat.i^ Pph, (27) 

where the summation index covers the different slices of atoms. With this Hamiltonian and the atomic decay and 
photon absorption loss mechanisms, the appropriate S,-, L,-, and N,- matrices are readily found, and the update 
formulae of Sec. ^ (or a slightly modified version thereof for the optically thick gas "2^) may be applied for the 
determination of the covariance matrix and the mean value vector for the Gaussian variables in y. The result of this 
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calculation is a 2n x 2n atomic covariance matrix, with only minor squeezing in each slice, as the quantum correlations 
are distributed over the entire sample. One readily obtains the noise properties of the total atomic spin components, 
but it is more interesting to find the smallest eigenvalue of the covariance matrix, corresponding to a specific spatial 
mode of the atoms which is maximally squeezed. This mode, indeed, is the one that couples most efficiently to the 
radiation, and it is hence this smallest eigenvalue that determines the precision with which one can estimate, e.g., the 
Larmor rotation rate of the collective SDin|2fi||. 



V. MAGNETOMETRY IN THE GAUSSIAN DESCRIPTION 



Precision atomic magnetometry relies on the measurement of the Larmor precession of a spin-polarized atomic 
sample in a magnetic field |2ll 122 . l2^ . From standard counting statistics arguments, one might expect the uncertainty 
in such measurements to decrease with the interaction time t and with the number of atoms A^at as 1/V-^ati- If, on 
the other hand, the monitoring of the atomic sample, necessary for the read-out of the estimate of the magnetic field, 
squeezes the atomic spin, the above limit may be surpassed. In a theoretical analysis '24'| it was suggested to estimate a 
scalar B field by a polarization rotation measurement of a far off-resonant light beam passing through a trapped cloud 
of spin-1/2 atoms. By quantum trajectory theorv psll co mbined with the classical theory of Kalman filters |2J) 113 j the 
uncertainty in the classical field strength was found '2 4^1 to decrease as l/(iVati^^^). This proposal was implemented 
experimentally, and indeed sub-shot-noise sensitivity was found 27]. In our analysis of the experiment 7, 2^, we 
advocated treating all variables, including the magnetic field, as quantum variables, and to assume a Gaussian 
probability distribution for the classical variable, so that the entire system can be described by the covariance matrix 
formulation. 

In the case of a scalar field directed along the y direction, the effective Hamiltonian of the system is given by 



(28) 



where /x^ = {t IK)[3\J \ ( J^) \/h is given by the magnetic moment /3, and where the B field causes a Larmor rotation of 
the atomic spin towards the z axis. Figure |3 shows the setup. It is the coupling of the B field to the spin-squeezed 
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FIG. 3: Setup for measuring the y-coordinate of a magnetic field. This is done by measuring the Farady rotation of a linearly 
polarized optical beam propagating through the atomic gas. 



variable Pat that makes an improved precision measurement of the magnetic field possible j29l|. 

The vector of variables in the case of a scalar magnetic field is y = (B, iatjPat, iph,Pph), and with the Hamiltonian 
of Eq. (EHll, the S'-matrix is found to be0,|2|| 
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(29) 



As By only causes rotation perpendicular to its direction, the variable iat oc Jy does not couple to {By,pat) and, 
hence, we only need to consider a 2 x 2 system with y — {By^p^.tY' ■ In the noise-less case, the system may now be 
propagated in time with the discrete update formula of Sec. |n] Alternatively we may consider the continuous limit 
and derive the differential equation for the covariance matrix A matrix of Eqs. (jS|)- (|12|l pertaining to y = (i?j,,pat)'^- 
The differential equation is on the matrix Ricatti form 8] 

A{t) = G - DA(i) - A(t)E - A{t)Y\{t), (30) 
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with G = 0, D = ( ° o), E = D^, and F = ([J ^2) where k"^ = k'^/t and /z = f-Lr/T. As may be checked by insertion, 
the solution to Eq. (1201) is = WU^i, where W = -DW + GU and U = FW + EU. The resuhing solution for 
the variance of the B field reads: 

Var(R(f^\ = Var(Bo)(K^t+l) /o-i \ 

\aL[^uyLj) iK4^2Var(So)f» + |KV^Var(Bo)t3+K2f+i ^'^'-^ 

—> 6 ^ 1 

The presence of noise'28j reduces the asymptotic decrease in the uncertainty with time from l/i'^ to 1/t. Figure 01 



1 - 




1 2 3 4 5 
t [ms] 

FIG. 4: Uncertainty of B field as a function of time. The value at f = 5 ms is ABy — 5.814 x 10"^ pT. We have chosen a 
segment duration r = 10~* s and corresponding field parameters = 0.0183 and /ir = 8.8 x 10"^. 

shows the decrease in the uncertainty of the B field with time in a calculation with physically realizable parameters. 

The new concept introduced in estimating the value of the classical B field is to treat the field itself as a quantum 
variable. Such an approach is not incompatible with the assumption that it is a classical parameter. We may imagine 
a canonically conjugate variable to B having an uncertainty much larger than required by Heisenberg's uncertainty 
relation and/or additional physical systems, entangled with the _B-variable, in which cases the S-distribution is 
indeed incoherent and "classical" . Also, one may argue that all classical variables are quantum mechanical variables 
for which a classical description suffices, and hence our theory provides the correct estimator according to the quantum 
theory of measurements: quantum mechanics dictates that the quantum state provides all the available knowledge 
about a system, and any estimator providing a tighter bound hence represents additional knowledge equivalent to a 
local hidden variable, and this is excluded by quantum theory. It is of course crucial that our measurement scheme 
corresponds to a quantum non-demolition (QND) measurement, i.e., we assume that there is not a free evolution of 
the B-field induced by its conjugate variable which may thus remain unspecified. It is also this QND property of the 
measurement scheme that implies a monotonic reduction of the uncertainty of B which is consistent with the classical 
parameter estimation (we can not unlearn what we have already learnt about B), unlike, e.g., the uncertainty of the 
atomic Xat variable which must increase when Var(pat) is reduced and when the atoms undergo spontaneous decay. 



VI. ENTANGLEMENT IN THE GAUSSIAN DESCRIPTION 

The theoretical proDOsal jSOl l3lj| and the subsequent experimental demonstration l32j that by using only coherent 
light, entanglement can be generated between distant free space atomic ensembles, has attracted much attention in 
the quantum information community. The primary reason being that entanglement in such macroscopic ensembles 
of atoms is robust and easy to make because of the effective and tunable atom-light coupling. Theoretical analyses 
of entanglement with continuous variables has been performed in the Heisenberg picture. 30| , and the state vector 
dynamics for a few tens of atoms has been considered by quantum trajectory methodspSl l34|. Also the Gaussian 
description has been successfully applied to describe the entanglement generation between two macroscopic atomic 
ensembles due to continuous probing of collective spin variables by optical Faraday rotation^ . 
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Generally, in the entanglement setup, the two gasses are polarized along opposite directions, say the positive and 
negative x axis. This means that the classical x components of the collective spin vectors are given by J^^i ~ 
^■^at,i/2 = Jx, and Jx.2 — —Jx, and the Gaussian description is applicable with the following vector of canonical 

quantum variables y = (xat,i,Pat,i, iat,2,Pat,2, iph,Pph)^ = (-^, --i^, ^)^. The 

Hamiltonian for either sample is given by Eq. I|21|l . To model the light-atom interaction, the light beam is divided 
into segments as discussed in Sec. Ill Dl The S"- matrices St,i and Si-,i for the two gasses are readily found from 
Heisenberg's equation of motion for the variables in y, and combined to St = St-,iSt-,2 for the full matrix. Additional 
evolution matrices may be defined that describe the rotation of the atomic variables of the samples and the effect of 
the homodyne detection jSj]. The update of the system then proceeds as outlined in Sec. ^ The theory incorporates 
the interaction between the atoms and the optical field, atomic decay, and the measurement induced transformation 
of the atomic state. The reduction of the full quantum state description to a simple Gaussian state fully represented 
by a set of mean values and a covariance matrix makes the system straightforward to deal with numerically, and 
analytical results can be obtained in several important cases. 

While the general problem of a measure for the entanglement between two mixed states remains unsolved, the 
entanglement between the atomic ensembles obtained by the continuous probing may be quantified by the Gaussian 
entanglement of formationls^ or the logarithmic negativity js^. The Gaussian description may also be used to identify 
the optimal performance of the entanglement scheme in the presence of atomic decay |^. 

A. Entanglement and vector magnetometry 

The possibility to entangle more atomic gasses was also considered in magnetometry [2II in connection with the 
problem of measuring two or three components of the B field with entangled gasses. In the case of two components, 
say By and Bz, the atomic sample is split in two and one gas is polarized along x and the other along —x. These 
polarizations assure that the two observables (Jj,j + Jy^) and (Jzj + Jza), and equivalently Xati — iat2 and pat,i +Pat2 
commute. (Note that a different sign convention for the fata variable was applied in our previous work|23.) The 
interaction between the magnetic fields and the two samples is described by the following effective Hamiltonian 

-ff^t.r = MrBj,(Xati - ^ats ) + MrB^feti +Pat2)- (32) 

This interaction causes changes in the atomic observables (pati — Pat2) and (xati +iat2) proportional with By and 
Bz, respectively. To probe these changes we introduce the effective light-atom interaction 

HLt.T = l^riPati - j5at2)Pphi + Kr(iati + iat2)^ph2, (33) 

where the appropriate relative sign between the atomic variables of the two gasses can be implemented by adjusting the 
sign on Kr after the probe beams have passed through the first gas 28J. The gasses are probed by the simultaneous 
action of the Hamiltonian from Eqs. H32l) - (|33|l . i?int,r — H^^^ ^ + Hl^-^^ ^. The vector of quantum variables is y = 

(Sz, Bj,, iati iPati 7 iat2 7Pat2 1 iphj jPphj , iph2 iPphj)"^- With this statc vcctor and the above Hamiltonian, the formalism 
of Sec. |n] can be directly applied and the final uncertainty of the B fields can indeed be lowered compared to the case 
with individual probe beams by letting the probe beams pass through both gasses and thereby entangling the two. 
An extension to full three dimensional vector magnetometry using three probing beams and six atomic samples can 
also be shown to have superior resolution in comparions with measurements on separable systems (28J. 

VII. EXTENSIONS OF THE THEORY 

In this section, we outline some topics which are subject to studies within the Gaussian description at the time of 
writing, and we discuss how to go beyond the Gaussian approximation. 

A. Non spin— 1/2 systems 

The theory presented in Sec. Illll explicitly used the representation of the collective angular momentum variable 
in terms of Pauli spin matrices = f This representation was crucial for the reduction of the Hamiltonian 

of Eq. 1)191) to the form of Eq. (|21|l which is expressed directly in terms of the canonical (x,p) Gaussian variables. 
For the more general problem of excited (\F' M')) and ground {\FM)) states with F > 1/2 and F' > 1/2, it is still 
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possible to obtain an approximate solution within the Gaussian description. This more complicated problem is of both 
fundamental interest and of practical importance since ongoing experiments based on the Faraday-rotation scheme 
are carried out with such real multilevel atoms [37ll8^ . 

Equation H17|l represents the interaction with atoms with arbitrary level structure, and to deal with such atomic 
samples within a Gaussian description, we suggest to introduce a second quantized formalism for the atoms in which 
bosonic atomic field operators ^m' create and destroy atoms with the given magnetic quantum number. Note 
that the bosonic character merely reflects the symmetry under permutations of the atoms: The theory works for both 
fermionic and bosonic atoms. We can then write the collective atomic operators in the Hamiltonian in terms of the 
atomic field operators, J^^ \FM)j{FM'\ = #j^^^I'^f/ to obtain 



Hint,r= J2 [(c+,M(^)a+a+ +C-,M(A)aLa-) *i/*M (34) 



At this point we make a mean field approximation, and we expand the field operators for the light fields 

d± a± + Sa±, (35) 

and the atom fields 

*A/ ^ $M + Si>M, (36) 

with c-numbers a± and $m, and "small" operators 6a± and S'^m- We insert Eqs. H35|) - I|3t)|) into Eq. H34|l and expand 
to second order in the operator terms. This procedure leads to a rather lengthy expression which is conveniently split 
into terms which are of zeroth, first and second order in the quantum fields. The classical fields are explicitly time- 
dependent and their dynamics is given by replacing all operators by their c-number parts in Heisenberg's equations of 
motion. Since we neglect operator terms above second order, the quantum part of the Hamiltonian is at most bilinear 
(with classical time-dependent coefficients). The linear terms cause mean drifts of the mean value of the operator com- 
ponents, which we can absorb in the c-number components. For the new quantum operators, we may then maintain 
{Sa±) — and = 0. In the resulting bilinear Hamiltonian involving the operator terms, we now make the tran- 

sition to the Gaussian state description by forming a vector of variables y — {x-f,P-f, ■ ■ ■ j 2;f,Pf, i;_ ,p_)"^, 

with x± = \[\(.5a± + 5a^^), p± = -i^{Sa± - Sal^), xm = \J^{5^ m + xm = -iy^((5^'j\/ - '5^'m), and 

adopting the formalism of Sec. ^ 

We note that the expansion of noise terms around classical mean values has been used as a standard tool in quantum 
optics, e.g., to deal with the optical Kerr-effect, and if only unitary dynamics and losses are considered, the present 
approach does not offer any new insights. It is important to remember, however, that we are also able to treat the 
dynamics conditioned on measurements on the system. Work is in progress, and we will report on results of this 
approach to the multilevel problem elsewhere. 



B. Quantum correlated light beams 



So far, we have treated the case of a coherent, monochromatic beam of light incident on the atomic samples. For 
high precision probing, atomic spin squeezing and entanglement, it has been proposed to use squeezed beams of light 
and twin beams, and this is indeed also possible within the Gaussian formalism. To model in a simple manner the 
coupling to squeezed light beams, one may simply alter the covariance matrix elements for the field operators in each 
beam segment prior to the interaction with the atomic sample, so that rather than the 2x2 identity matrix with 
equal variances of the two field quadratures in Ea. (|ll|l . we assume the form 

B = diag(l/r,r), (37) 

where r is the squeezing parameter. Carrying out the calculations as described in the previous sections, we observej3, 
|2^ . that the magnetometer resolution is improved by this parameter. 

As pointed out in our analysis^, however, a squeezed beam of light carries correlations between the the field 
operators evaluated at different times. This implies, that if one observes the beam for only a very short time, one 
will not be able to detect the squeezing. It is well-known from the standard quantum optical analysis of the optical 
parametric oscillator (OPO) (in the Heisenberg picture in the frequency domain) that there is a certain frequency 
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band-width, F, of squeezing associated with the field decay rate of the cavity. Only the accumulated signal over 
sufficiently long times {t > 1/r) will be able to extract the frequency components for which the field is squeezed. One 
might think, that one should therefore consider a longer string of light segments with appropriate covariance matrix 
elements, and carry out the update on all segments (and the atoms) every time a single optical segment is detected. 
In fact, there is an easier approach. The squeezed beam is produced by continuous leakage of the field inside the 
OPO cavity, and the temporal correlations are due to this joint source of the radiation. The most economical way 
to describe the interaction of atoms with a squeezed beam is therefore to incorporate the single-mode field inside 
the cavity in the Gaussian state formalism, and to consider again only one optical beam segment at a time, from its 
creation out of the cavity, interaction with the atoms, and final detection, and hence causing an update of the joint 
atom and cavity covariance matrix. 

We have implemented such a modeljs^ and verified that it reproduces the known noise properties for the signal 
integrated over both short and long times. In addition, we have applied the model to magnetometry, and observed 
that for segments shorter than the inverse band- width of squeezed light, the resolution is not improved with the use 
of squeezed light, but after many segments and a long total interaction ti me, the results asymptotically approach the 
factor 1 /r improvement of the simple model as one might have expected • 

In closing this section, we note that the Gaussian description is not restricted to the examples and extensions 
discussed above. Extra physical systems can be included straightforwardly by adding appropriate rows and columns 
to the covariance matrix. In this way, one may, e.g., describe the effects of imperfect detectors and filters, and loss in 
optical fibres. 



C. Beyond the Gaussian approximation 

Looking back on the development of the theory in this work, we note that the interaction and the dissipation can 
be treated exactly without recourse to a Gaussian ansatz for the quantum state. In fact, Eq.(|Sll is a general update 
formula for the Wigner function under an arbitrary measurement, but in the general case this expression may be 
difficult to evaluate, and in particular to use as input in the next step of the continuous probing on the system. In a 
recent experiment |40j. a Gaussian squeezed state was mixed with the vacuum field at a beam splitter to produce an 
entangled two-mode field state part of which was sent to an avalanche photodiode and part of which was monitored 
by homodyne detection. The state of the second component conditioned on a photodiode countin g ev ent is a non- 
Gaussian state, as verified by a double-peaked homodyne detection signal. The process was modeled|41| by assuming 
that prior to the registration of a single photon, the state of the field is described by the Gaussian Wigner function 
VVcausslT, <5) with 7 and 5 denoting pairs of real variables of the two modes. The detection of a single photon, |1), 
corresponds to apphcation of the Wigner function W,neas{S) = W|i)(i|(^) = f exp(— 2|(5p)(— 1 + 4|(5p) in Eq. |(SJl,and 
the Wigner function for the second beam conditioned on this state, Wc{j) = J W(7, S) yV|i)(i| (S) (PS is readily shown 
not to be a Gaussian. Another example where one 'jumps' out of the Gaussian states is in a recent proposal 42J where 
a setup of two beam splitters with carefully chosen properties and photodetectors allows one to produce with high 
fidelity a single-photon state from Gaussian squeezed vacuum input beams. We believe that a combination of the 
theory of Gaussian state updates and inclusion via Eq. (|SJ) of one or a few non-Gaussian preserving measurements 
may be a useful approach to these problems. 

As it is necessary to leave the Gaussian states to perform some quantum information tasks such as distillation 
of entangled states^, 0, it is in general important to have tools to handle the interface between Gaussian and 
non-Gaussian states. 



VIII. CONCLUSIONS AND OUTLOOK 



We have presented a general formalism to treat the dynamics of 2n canonical variables within a Gaussian description. 
For a probing light beam, the concept of "segment quantization" where the beam is quantized in small fragments of 
duration r and length L — ct allows an efficient description of not only the evolution of the system subject to the 
Hamiltonian, but also to the measurement process through simple update formulae for the mean value vector and the 
covariance matrix, which fully characterize the Gaussian state. Of particular current interest in the field of quantum 
information and quantum communication [loll3Cll3^ . and in precision magnetometrvp^ . is the off-resonant probing of 
ensembles of atoms leading to a dispersive Faraday effect. This interaction is bilinear in the effective canonical variables 
of the system^nd a description within the Gaussian framework of, e.g. spin squeezing j205, magnetometrv,?,. .28iJ . and 
entanglement Q is straightforward. 

In Sec. lVllI we have outlined some possible extensions to the Gaussian description. In the future it will be interesting 
to develop further theory for continuous variable system which leave the Gaussian description, either because of the 
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interactions involved, because of the measurement schemes, or because of coupling of a small discrete system to 
the collective continuous degrees of freedom (examples: single photons can be stored and emitted on demand by 
macroscopic atomic samples, trapped ions can be entangled by continuous probing with classical laser fields). Such 
approaches hold the potential to form Schrodinger Cat states, which may have favorable properties in high precision 
detection, and they may be used to implement distillation and purification protocols from quantum information theory, 
which are known not to work for Gaussian states and operations. The theoretical task is to identify processes that can 
be implemented experimentally and which break the Gaussian character, and to establish a theoretical description of 
the resulting states, which will invariably be much more complicated to deal with than the Gaussian states. 
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